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1.  Introduction 


Several  algorithms  have  been  described  in  the  literature  for  the  numerical  representation  of  the 
inviscid  flux  in  the  Navier-Stokes  equations.  The  more  recent  of  these  schemes  have  their  un¬ 
derlying  basis  in  either  mathematical  (for  example  TVD  schemes  [1])  or  physical  (for  example 
flux-splitting  [2,  3])  analyses.  All  of  these  modern  methods  employ  upwinding  in  one  form  or 
another  to  obtain  algorithms  possessing  better  dissipation  characteristics,  higher  stability  bounds 
and  increased  numerical  efficiency  [4,  5].  This  inherently  dissipative  nature  of  upwinded  schemes 
is  beneficial  at  or  near  shock  waves  where  the  inviscid  fluxes  dictate  finite  jumps  in  flow  variables. 
However,  typical  viscous  supersonic  flows  about  bodies  display  steep  gradients  not  only  in  shock 
regions  but  also  in  the  boundary  layers.  The  dissipative  characteristics  of  upwind  schemes  in  these 
latter  regions  is  a  drawback  since  the  numerical  dissipation  may  be  comparable  to  or  overwhelm 
the  viscous  terms.  This  can  often  have  detrimental  consequences  on  the  prediction  of  heat  transfer 
rates  [6]. 

Following  the  approach  of  several  authors  (7,  8]  upwind  schemes  may  be  classified  according 
to  the  Riemann  solver  used  to  evaluate  the  fluxes  at  the  cell  faces  and  the  mechanism  through 
which  higher  order  spatial  accuracy  is  obtained  —  non-MUSCL  or  MIJSCL.  Some  examples  of 
schemes  obtaining  higher  order  accuracy  with  non-MUSCL  approaches  are  the  “Symmetric”  TVD 
scheme  of  Yee  [9,  10]  and  the  “Upwind”  TVD  scheme  of  llartcn  and  Yee  [10,  11]  where  the  quotes 
indicate  the  fact  these  words  do  not  retain  the  traditional  centered  or  upwind  meanings  because 
of  the  presence  of  flux  functions  and  slope  limiters. 

The  focus  of  the  present  research  is  on  schemes  that  obtain  higher  order  accuracy  through  the 
MUSCL  approach  of  van  Leer  [12].  Specifically,  we  study  the  performance  of  the  flux  vector  s|)lit 
methods  of  MacCormack  and  Candler  [6]  (to  be  abbreviated  henceforth  as  the  MC  method)  and 
van  Leer  [3]  (abbreviated  vL)  and  the  flux  difference  split  method  of  Roe  [13].  The  distinction 
betweeti  flux- recto/- and  flux  (//jffc /■(  ///•/  sjilitting  may  be  based  upon  the  model  used  to  ilerive  the 
basic  interaction  between  neighboring  cells.  When  the  interaction  between  neighboring  cells  is 
mod/'led  with  finite- ami)litude  waves  (the  Riemann  approach),  the  scheme  is  called  flux-difference 
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splitting;  when  accomplished  through  the  mixing  of  pseudo-particles  (the  Boltzmann  approach) 
the  method  is  called  flux- vector  splitting  [3]. 

The  original  Steger-Warming  flux-vector  split  method  [2]  (often  called  the  beam  scheme) 
exploits  the  homogeniety  of  the  inviscid  flux  vector  to  split  the  flux  into  positive  and  negative 
components  depending  upon  the  sign  of  the  local  eigenvalues.  Since  the  mass  flux  for  this  scheme 
is  not  continuous  where  the  eigenvalues  change  sign,  “glitches”  are  often  observed  at  sonic  and 
stagnation  points.  This  method  was  proven  to  be  excessively  dissipative  in  the  boundary  layers  by 
MacCormack  and  Candler  [6]  who  suggest  relatively  simple  changes  (outlined  later)  to  eliminate 
numerical  diffusion;  these  modifications  give  rise  to  the  MC  method.  Reported  applications  of  the 
original  Steger-Warming  algorithm  include  inviscid  flows  past  airfoils  [2]  and  cylinders  [14].  The 
MC  scheme  has  previously  been  applied  for  flat  plate  boundary  layer  flows,  flows  past  compression 
ramps,  blunt  body  flows  [6],  viscous  real  gas  flows  past  sphere-cones  [15],  Type  III'^  and  Type  IV 
viscous  shock-shock  interactions  [16]. 

Unlike  Steger  Warming  splitting,  van  Leer’s  formula  exhibits  flux  continuity  across  eigenvalue 
sign  changes.  Some  improvement  is  obtained  in  shock  capture  capability  over  the  Steger  Warm¬ 
ing  method  [17]  .  Applications  with  van  Leer’s  scheme  include  inviscid  subsonic  and  transonic 
flows  over  airfoils  [18],  viscous  shock-induced  separated  flows  [19],  flows  over  delta  wings  [20] 
and  recently  a  Type  IIl'^  interaction  at  Mach  8  [21]  with  and  without  a  turbulence  model.  A 
comparison  of  the  van  Leer  and  the  original  Steger-Warming  scheme  for  some  Euler  flows  may  be 
found  in  Anderson  et  al.  [17]. 

Roe’s  flux-difference  split  method  [13]  is  based  upon  accurate  prediction  of  wave  interactions 
between  interfaces  through  an  approximate  (linearized)  equation.  This  approach  has  been  demon¬ 
strated  to  reduce  the  numerical  dissipation  [22]  and  has  subsequently  been  extended  to  reduce 
grid-dependence  [23]  common  to  methods  based  on  1-D  analysis.  One  drawback  of  this  algorithm 
is  the  possible  violation  of  the  entropy  condition,  necessitating  the  employment  of  appropriate  en¬ 
tropy  correction  factors.  One  approach  is  described  in  considerable  detail  in  Chapter  2.  Roe’s  flux 
difference  split  method  has  been  applied  in  recent  years  for  viscous  and  inviscid  conical  flow  [13], 
nozzle  flow  and  shock  reflection  [24]. 
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The  popularity  of  the  above  three  flux  split  methods  in  the  computation  of  high  speed  flows 
is  likely  to  increase  especially  with  their  extensions  to  include  high  temperature  effects  by  Liu 
and  Vinokur  [25],  Grossman  and  Cinnella  [26,  27]  and  Shuen  et  al.  [28].  A  comparison  of  the 
relative  merits  of  each  of  these  schemes  has  been  performed  either  on  inviscid  flows  [17,  8]  or 
on  simple  viscous  flows  [22,  29],  The  former  type  of  comparison  focuses  exclusively  upon  shock 
capture  capability  while  the  latter  is  restricted  to  relatively  simple  flows  (e.g.,  flat  plate  boundary 
layers).  An  overall  assessment  of  these  methods  based  upon  the  studies  available  in  the  literature 
is  further  complicated  by  the  ambiguity  in  specific  implementation. 

We  perform  a  detailed  analysis,  with  grid  resolution  study,  of  the  prediction  capabilities  of 
these  three  algorithms  on  two  classic  2-D  viscous  laminar  problems  with  specific  emphasis  in 
the  prediction  of  heat  transfer  rates.  The  implementation  of  each  algorithm  is  determined  for 
our  purposes  to  be  the  one  most  commonly  utilized  as  observed  in  the  literature.  The  choice  of 
problems  and  flow  parameters  is  dictated  by  the  need  to  simulate  features  generic  to  a  wide  range 
of  complex  flows  for  which  reliable  experimental  data  in  the  form  of  surface  pressure  and  heat 
transfer  are  available. 

•  Mach  16  Blunt  body  flow:  This  is  a  classic  problem  typified  by  a  simple  geometry  and 
flow  features  and  the  availability  of  a  relatively  accurate  theoretical  value  of  stagnation 
heat  transfer  [30].  Despite  the  relative  simplicity  of  the  flowfield,  several  algorithms  have 
displayed  difficulties  in  heat  transfer  prediction.  These  difficulties  generally  manifest  them¬ 
selves  as  the  so-called  “carbuncle”  or  “bulge”  phenomenon  [31,  32,  33]  and  may  be  the  result 
of  the  singularity  of  the  eigenvalues  across  the  entire  length  of  the  line  of  symmetry. 

•  Mach  14  flow  past  a  compression  ramp:  This  type  of  flow  exists  generically  in  several 
practical  applications  (e.g.,  inlets)  and  is  typical  of  shock-wave  boundary  layer  interactions 
(Fig.  I ).  The  parameters  chosen  duplicate  one  of  the  experiments  of  Holden  and  Moselle  [34] 
corresponding  to  a  24°  ramp  at  Mach  14.1.  The  resultant  flowfield  is  known  to  exhibit  a 
large  region  of  recirculation  [3.5].  Although  some  evidence  exists  to  suggest  that  3-D  effects 
may  be  important  [36],  no  effort  is  made  to  resolve  this  issue  in  the  present  work. 
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In  each  instance,  a  sequence  of  grids  is  employed  to  obtain  grid  resolution  studies.  The  finest 
grids  utilized  are  presented  in  Fig.  2.  Further  quantitative  details  are  presented  and  discussed 
with  the  results. 
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14.1 


(b)  Corner  flow 


Figure  2:  Grids  employed 
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2.  Theoretical  Model 


The  2-D  Navier  Stokes  equations  in  strong  conservation  form  are  solved  in  the  transformed  (^,  t?) 


coordinates; 


^  ^  ^  -  n 

^  ^  dr} 


where  U  is  the  solution  vector  of  conserved  variables  {{p,  pu,  pv ,  pe)  j  J)  with  p  denoting  density,  u 
and  V  the  Cartesian  components  of  velocity,  e  the  total  energy  per  unit  mass  (=  e,  +  0.5(u^  +  v^)) 
and  ei  the  internal  energy  per  unit  mass  (=  C„T)  where  Cy  is  the  specific  heat  at  constant 
volume  and  T  is  the  static  temperature.  The  flux  vectors,  F  and  G,  include  the  full  viscous  and 
inviscid  terms,  and  T](x,y)  are  the  transformed  variables.  For  the  configurations  under 

consideration,  the  ^  coordinate  lines  are  in  the  generally  circumferential  or  streamwise  direction 
while  the  t]  lines  are  radial  or  body-normal.  The  density,  static  pressure  p  and  temperature 
are  related  through  the  equation  of  state  p  =  pRT  where  R  is  the  gas  constant.  The  molecular 
viscosity  p  is  approximated  by  Sutherland’s  law  and  the  molecular  Prandtl  number  Pr  is  assumed 
to  be  0.73  (air). 

Equation  2.1  may  be  interpreted  as  describing  the  balance  of  mass,  momentum  and  energy 
over  an  arbitrary  control  volume  [18].  Numerical  time  integration  is  achieved  with  a  residual 
driven  line  Gauss-Seidel  relaxation  scheme  as  described  in  the  works  of  MacCormack  [37]  and 
Candler  [38].  With  first-order  backward  Euler  time  discretization  and  linearization  of  the  fluxes 
in  time,  the  discretized  equation  may  be  written  in  the  form: 

{NUMER1CS}6U  =  PHYSICS  (2.2) 


where  PHYSICS  represents  the  residual  and  NUMERICS  contains  the  driving  terms  and  6U  rep¬ 
resents  the  change  in  the  solution  vector  at  each  time  step.  The  full  Navier-Stokes  equations  are 
utilized  in  computing  the  residual  with  the  appropriate  upwind  model  for  the  inviscid  terms  and 
centered  evaluation  of  the  viscous  terms.  One  advantage  of  the  above  methodology  is  that  approx¬ 
imations  may  be  utilized  in  the  NU M ERICS  portion  of  the  code  without  affecting  the  accuracy 
of  the  converged  result.  In  fact,  Liou  and  van  Leer  [24]  point  out  that  even  if  the  PHYSICS 
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term  is  evaluated  with  other  methods  (such  as  van  Leer’s  splitting  or  Roe’s  upwinding),  the  use 
of  Steger  Warming  splitting  in  evaluating  the  NU M ERICS  portion  leads  to  robust  codes  for  the 
Newton  linearization  procedure.  The  implicit  operator  in  this  research  therefore  utilizes  the  Ste¬ 
ger  Warming  Jacobians  to  obtain  strong  diagonal  dominance.  Further,  for  simplicity  the  viscous 
Jacobians  in  the  driving  NU M ERICS  portion  are  computed  with  the  thin  layer  approximation. 

When  solved  with  the  Gauss-Seidel  line  relaxation  procedure,  Eqn.  2.2  represents  a  block 
tridiagonal  system.  For  transonic  and  supersonic  flows,  line  relaxation  methods  are  superior  to 
approximate  factorization  methods  in  convergence  rate  which  amply  compensates  for  the  larger 
computation  required  per  iteration  [39].  The  line  Gauss-Seidel  algorithm  is  also  unconditionally 
stable  in  the  linear  analysis  and  is  known  to  be  relatively  insensitive  to  the  choice  of  time  increment 
per  iteration. 

For  most  computations,  the  automatic  CFL  number  adjustment  algorithm  of  MacCormack  [40] 
is  utilized.  In  this,  the  allowable  time  step  is  automatically  permitted  to  double  every  nq  iterations 
from  an  initial  value  (typically  small  ~  0.01).  After  the  explicit  part  of  each  iteration  (the  first 
stage  if  a  two-stage  method  is  used),  the  changes  in  the  solution  are  adjusted  so  that  the  maximum 
relative  change  in  density  at  any  point  due  to  the  present  iteration  is  always  less  than  a  fixed 
value  {Spniax  say)  and  the  maximum  relative  decrease  in  temperature  is  also  less  than  another 
fixed  value  (^T„ii„  say).  This  adjustment  involves  the  division  of  the  residual  at  each  point 
by  a  fixed  number  determined  such  that  the  above  conditions  are  satisfied.  In  effect,  this  cuts 
the  explicit  allowable  timestep  by  cutting  the  CFL.  Although  MacCormack  obtained  arbitrary 
increases  in  the  CFL  number,  for  some  of  the  current  computations,  the  residuals  as  well  as  CFL 
numbers  asymptote  to  limit  cycles.  In  such  instances,  the  maximum  CFL  is  restricted  to  such 
a  value  that  the  CFL  number  remains  stable  and  the  solution  residual  converges.  This  typically 
eliminates  the  high  frequency  oscillations  in  the  residual  cycles.  Also,  although  MacCormack 
used  maximum  allowable  relative  change  values  {\Sp/p\max  or  -ST/T\jnax)  of  0..5,  in  the  present 
instance,  especially  for  finer  grids,  values  of  0.01  to  0.05  were  used. 

A  brief  description  of  each  algorithm  is  presented  with  reference  to  the  flux  evaluation  ((j)at 
a  j  -f  I  surface  in  the  present  cell  centered  finite  volume  formulation.  Denoting  the  derivatives 
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of  the  transformed  coordinates  with  respect  to  the  physical  coordinates  as  etc.,  the  flux  at 

this  surface  in  generalized  coordinates  is: 

G  =  J-^  {vxF  +  rjyG)  (2.3) 

where  J  is  the  Jacobian,  F  and  G  are  the  fluxes  expressed  in  Cartesian  coordinates  and  include 
viscous  and  inviscid  components: 

F  =  F  +  and  G  =  G  +  (2.4) 

with  F  and  G  representing  Cartesian  inviscid  fluxes  and  F-u  and  G„  the  viscous  fluxes  respectively. 
The  viscous  fluxes  are  evaluated  in  a  centered  fashion  and  will  not  be  discussed  extensively.  The 
vectors  grad(^)/J  and  grad(7/)/J  represent  the  directed  areas  of  cell  interfaces  in  the  rj  and 
C  directions  and  J  is  the  inverse  ceU  volume.  Since  the  extension  to  general  coordinates  is 
straightforward  as  described  by  Anderson  et  al.  [17],  the  following  discussion  focuses  only  on 
the  inviscid  Cartesian  flux  vector  G. 

At  each  interface,  the  state  of  the  flow  is  described  by  two  vectors  of  conserved  variables, 
and  on  either  side  of  the  interface.  These  vectors  are  obtained  from  the  known  values  at  the 
cell  centers  with  the  MUSCL  approach  in  conjunction  with  a  limiter  as  described  below. 

2.1  MacCormack  and  Candler  Flux- Vector  Split  Algorithm 

In  flux  vector  splitting,  the  inviscid  flux  is  split  into  positive  and  negative  components  for  appro¬ 
priate  upwind  differencing: 

=  +  (2.5) 

We  describe  first  the  original  Steger  Warming  method  and  then  the  MacCormack  and  Candler 
scheme.  In  the  Steger  Warming  approach  [2]  (abbreviated  OSW),  since  the  inviscid  fluxes  are 
first  order  homogeneous  and  hyperbolic,  the  flux  Jacobian  B  {B  =  dGfdU,  G  =  BU)  may  be 
written  as: 

B  =  Q-^AQ  =  Q-‘(A+  -h  A-)Q  =  5+  -b  B"  (2.6) 
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where  A  is  a  diagonal  matrix  consisting  of  the  eigenvalues  of  (A  =  diag{t;,t)  +  c,v-  c})  and 
A"*"  and  A“  denote  the  splitting  of  the  eigenvalues  into  positive  and  negative  components.  For 
convenience,  the  matrix  Q  may  be  written  as  [38]: 

Q  =  CS  (2.7) 

where  S  =  V  is  the  vector  of  nonconserved  variables  {p,  u,n,p}  and  C  diagonalizes  the  flux 
vector  G  written  in  terms  of  V.  At  a  face  j  +  1/2  therefore,  the  inviscid  original  Steger  Warming 
flux  (Gosw)  may  be  written  as: 

Oosw.,^i  =  BfU‘-  +  B-^,V'<  (2.8) 

For  second  order  accuracy,  the  vectors  of  conserved  variables,  f/^and  U^,  are  obtained  at  the 
cell  interfaces  by  extrapolation  to  the  cell  surface  with  the  MUSCL  approach  of  van  Leer  [12]  in 
conjunction  with  the  minmod  limiter.  Denoting  the  vector  of  extrapolated  quantities  as  W: 

=  Wj„  -  (2-9) 

=  +  (210) 

where: 

=  minmod(A^^i,  Aj_i)  (2.11) 

and  A  ,  1  =  .  Several  choices  of  W  are  possible,  e.g.,  the  vector  of  conserved  variables 

2 

(W  =  U)  OT  characteristic  variables.  For  the  present  method  we  choose  to  extrapolate  the 
primitive  variables  {W  =  {p,u,v,p))  and  U  is  extracted  from  the  extrapolated  W.  With  the 
addition  of  the  minmod  limiter,  the  algorithm  reverts  to  first  order  accuracy  at  shocks  in  order 
to  preserve  monotonicity  within  the  solution. 

The  formula  of  Eqn.  2.8  exhibits  two  major  difficulties.  First,  it  may  introduce  discontinuities 
in  the  solution  at  sonic  and  stagnation  points  where  the  eigenvalues  change  sign  and  some  form 
of  eigenvalue  smoothing  may  be  necessary.  A  more  serious  flaw  in  Eqn.  2.8  is  that  it  introduces 
excessive  numerical  damping  in  the  boundary  layers.  This  damping  significantly  deterioriates 
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the  prediction  accuracy  of  surface  quantities  of  engineering  interest.  With  direct  algebraic  ma¬ 
nipulation  of  the  normal  flux  near  a  surface  under  boundary  layer  conditions,  MacCormack  and 
Candler  [6]  proved  that  the  Steger  Warming  procedure  introduces  artiflcial  tangential  momentum 
exchange  between  adjacent  cells  in  the  boundary  layer  solely  due  to  the  splitting  of  the  inviscid 
fluxes.  This  diffusive  term  is  proportional  to  the  quantity: 

~  “f,j)  (2.12) 

where  c  is  the  speed  of  sound  and  7  is  the  ratio  of  specific  heats.  Since  it  is  of  order  Ay,  it  can 
obtain  unacceptably  high  values  in  the  boundary  layer.  In  addition,  there  is  also  a  large  numerical 
exchange  of  kinetic  energy  of  the  order: 

~  ^i,j)  (2.13) 

between  adjacent  points  in  the  boundary  layer  where  a  is  the  kinetic  energy.  They  recommended 
that  both  and  B~  ,  be  evaluated  at  the  same  point.  Two  choices  are  obvious:  1)  evaluate 

J'T  2  J'r  2 

5  at  j  -f  5  by  averaging  the  conserved  variables  between  points  j  and  j  -I-  1,  or  2)  evaluate  B 
at  j  and  j  -f  1  respectively  in  successive  iterations.  No  significant  difference  is  observed  between 
results  obtained  with  either  approach.  In  the  present  calculations,  the  first  approach  is  utilized 
exclusively.  The  inviscid  MacCormack  and  Candler  flux  {Gmc)  reads  therefore: 


'MC,]  + 


.  =  Bl,U^  +  B-.U^ 


=  5-‘C-*|A+|C5]  -I-  [5-*C-‘|A-|C'5]  U 


(2.14) 


where  once  again,  and  are  evaluated  with  the  MUSCL  approximation  (Eqns  2.9,  2.10 
and  2.11).  This  modification,  however,  introduces  an  artificial  pressure  gradient.  This  error 
comprised  the  terms: 


(2.15) 


and  may,  therefore,  be  expected  to  dominate  only  in  the  close  vicinity  of  the  boundary  where 
high  velocity  gradients  exist.  MacCormack  and  Candler  [6]  treated  this  problem  at  the  expense 
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of  reintroducing  some  of  the  diffusiveness  of  the  original  Steger  Warming  scheme  by  further 
modifying  the  component  5  for  the  matrix  Q  in  Eqn.  2.7  such  that  its  last  row  is  evaluated  as 
in  the  original  formulation.  Subsequently,  MacCormack'  suggested  that  since  this  error  arises 
due  to  the  multiplication  of  the  last  row  of  the  matrix  5  and  the  column  vector  U ,  this  product 
be  explicitly  replaced  with  the  extrapolated  pressure  at  the  interface.  This  latter  approach  is 
followed  in  the  present  computations. 

The  MC  method  was  developed  for  application  only  in  the  boundary  layers  and  in  fact  leads 
to  instability  when  applied  in  regions  of  discontinuities  such  as  shock  waves.  As  a  result,  it  is 
necessary  to  revert  to  the  original  Steger  Warming  scheme  near  shocks.  This  is  achieved  in  a 
smooth  and  systematic  manner  by  defining  the  parameter 


Xi  = 
PR  = 


1 

1  +  PRx  PR 
\Pi,j+\  -  P«.j  l 


and  defining  the  flux  at  the  j  +  1/2  face  as: 


(2.16) 

(2.17) 


Gj+L  =  +  (1  -  Xi)Gosw,/^^^ 


(2.18) 


2.2  van  Leer’s  Flux- Vector  Split  Algorithm 

The  functional  form  of  the  van  Leer  scheme  [3]  is  similar  to  the  Steger  Warming  algorithm.  For 
supersonic  flow: 

G+  =  G;  G~  =  0  for  My  >  1  (2.19) 

G+  =  0;  G-  =G  for  My  <  -  I  (2.20) 

where  My  is  the  local  Mach  number  normal  to  the  face  j  +  For  subsonic  flow,  G  was  revised  by 
van  Leer  to  avoid  the  discontinuity  exhibited  across  sonic  lines  by  the  Steger  Warming  algorithm. 

‘Private  Communication,  June,  1990 
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2.3  Roe’s  Flux- Difference  Split  Algorithm 

The  formula  for  this  scheme  reads  [13]: 

(t'Ai  -  f'i.i)  (2-24) 

where  ( ' )  indicates  evaluation  at  the  Roe  averaged  state  between  and  (7^[13].  This  scheme 
extends  to  second  order  accuracy  through  the  MUSCL  approach  (Eqns.  2.9  and  2.10)  with  W  = 
{p,u,v,p}.  In  contrast  to  the  above  flux  vector  split  algorithms,  Roe’s  scheme  may  violate  the 
entropy  condition  when  the  eigenvalues  at  the  Roe  averaged  state  vanish.  Following  Harten  [41] 
the  eigenvalues  |A|  of  |v4l  are  modified  when  they  are  below  some  small  threshold  6: 

|A|  =  l-^l  <  ^  (2.25) 

'.vher^'  th?  value  of  S  is  taken  as  the  following  [10]: 

S  =  6J-^  ju.V^I  +  |«.Vr?|  +  ^  (|V^|  +  |V^|)  (2.20) 
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The  value  of  S  utilized  typically  in  the  range  0.05  to  0.1.  For  the  flow  past  the  compression  ramp, 
only  the  two  eigenvalues  of  the  form  u  +  c  and  u  -  c  (corresponding  to  the  genuinely  nonlinear 
eigenvectors)  are  cutoff  in  each  direction.  For  blunt  body  flows,  however,  it  is  necessary  also  to 
apply  the  cutoff  algorithm  to  the  eigenvalue  of  the  form  u  (corresponding  to  the  linearly  degenerate 
eigenvector)  in  the  streamwise  (body- tangential  or  ^)  direction  for  stability  purposes  [42].  Further, 
the  value  of  the  cutoff  parameter,  S  also  needs  to  be  increased  to  prevent  the  development  of 
anomalous  solutions  as  discussed  later.  This  can  however  lead  to  excessive  dissipation.  To  prevent 
this  for  both  the  blunt  body  and  the  corner  flow,  an  anisotropic  cutoff  formula  described  by 
Muller  [43]  is  employed  in  the  streamwise  ^  direction: 


S  = 


(2.27) 


where  =  |u.VA:|  Further  discussion  on  the  choice  of  the  parameter  6  is  provided  with 

the  results. 
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3.  Boundary  Conditions  and  Numerical  Details 


For  each  flow,  the  boundaries  may  be  categorized  into  one  of  the  following: 


•  “Inflow/Farfield  boundary:”  The  flow  vector  {p,  pu,  pv,  pe}  is  specified  corresponding  to  the 
known  values. 

•  Solid  boundary:  The  velocity  vector  and  the  normal  pressure  gradient  are  assumed  zero  and 
a  fixed  surface  temperature  is  specified  i.e.: 

dp 


pu  =  0;  T  =  Tu,; 


dn 


=  0 


(3.1) 


where  v  is  the  velocity  vector  and  the  subscript  w  refers  to  wall  conditions.  The  values  of 
Tu,  are  provided  with  the  results. 

•  Outflow  boundaries:  The  flow  at  these  boundaries  is  assumed  predominantly  supersonic. 
The  zero  gradient  extrapolation  condition  (d/d^  =  0)  is  applied. 

The  boundary  conditions  for  the  implicit  portion  of  the  algorithm  are  described  in  Gaitonde  and 
Shang  [16]. 

Convergence  is  determined  for  all  computations  by  monitoring  several  quantities.  The  global 
residual,  defined  as: 

1 


1IG./2.II  = 


IL  JL  4 

EEE 


Rk 


(3.2) 


is  one  measure.  In  Eqn.  3.2,  k  denotes  the  equation  ( 1  =continuity,  2,3=momentum  and 
4=energy)  on  an  IL  x  .1 L  computational  mesh  and  Rk  is  the  residual: 

dUk  dPk  dGk 


Rk  = 


(3.3) 


dt  d(,  dr] 

For  plotting  purposes,  these  values  are  normalized  by  the  value  of  ||G'./?.||  obtained  after  the 
first  iteration.  With  this  criterion,  convergence  is  eussumed  after  ||G./Z.||  drops  8  or  more  orders 
of  magnitude.  After  this,  the  surface  pressure  and  heat  transfer  are  monitored  over  about  one 
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characteristic  time  (Tc  =  LjUoo  where  Uao  is  the  freestream  velocity  and  L  is  a  macroscopic  length 
scale,  the  diameter  of  the  cylinder  for  the  blunt  body  flow  or  the  length  from  the  leading  edge 
to  the  corner  for  the  corner  flow).  In  addition,  the  integrated  root  mean  square  (RMS)  pressure 
and  heat  transfer  values  over  the  entire  surface  are  also  monitored.  Both  must  remain  constant 
at  convergence.  With  p  denoting  pressure  the  RMS  surface  pressure  may  be  written  as: 


RMSP  = 


1 


IL  , 

I  Pi,}=aurfa 

.=  l  ^ 


With  Q  denoting  surface  heat  transfer,  the  RMS  Surface  Heat  transfer  is: 

1 


IL 


RMSQ  =  — ,  5];(Q..j=.„r/ace)' 


(3.4) 


(3.5) 
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4.  Blunt  Body  Flow 


4.1  Flow  Parameters  and  Grid  Details 

The  parameters  for  this  flow  are; 

Moo  =  16.34 
Cylinder  Radius  =  1.5  in 


Poo 

Reynolds  Number 
^wall 


93.93R 
0.01203psia 
1.2  X  lO^per  foot 
530R 


This  case  represents  a  low-enthalpy  flow  for  which  previous  computations  validate  the  perfect 
gas  approximation  [44]).  Three  grids  (denoted  1,  2  and  3  respectively)  are  developed  with  the 
characteristics  described  in  Table  4.1.  In  this  table,  the  number  of  flow  points  listed  in  each 
direction,  IL  X  JL,  do  not  include  shadow  cells  employed  on  the  boundary.  Guidelines  for  grid 
resolution  are  taken  from  the  work  of  Klopfer  and  Yee  [45]  who  recommend  a  surface  cell  Reynolds 
number  (Rcc)  of  roughly  3  for  heat  transfer  calculations  (fixed- wall  temperature)  and  of  the  order 
of  10  for  adiabatic  wall  conditions  for  their  TVD  scheme.  They  also  indicate  that  heat  transfer 
rates  are  not  particularly  sensitive  to  the  circumferential  spacing.  The  most  refined  grid  for 
each  computation  satisfies  these  criteria  and  is  generated  with  a  combination  of  exponentially 
stretched  and  uniform  spacings  (Fig.  2(a)).  From  this  complete  grid,  subgrids  1  and  2  are  extracted 
systematically  to  provide  grid  refinement  studies.  Note  that  the  surface  cell  Reynolds  numbers 
presented  in  Table  4.1  are  based  on  the  half-cell  heights  of  the  first  cell  on  the  body. 
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Table  4.1:  Blunt  body  grid  details 


Grid 

IL 

X 

JL 

^^dmin 

R^dmax 

Re  1 

^^Imax 

A0|a. 

m 

1 

19 

X 

17 

19.4 

125.3 

51.5 

1.10 

10.5 

6.6 

1.98 

2 

39 

X 

35 

8.1 

51.6 

20.6 

0.51 

5.3 

3.4 

1.41 

3 

79 

X 

71 

IQII 

24.2 

9.6 

0.25 

2.5 

1.7 

1.18 

Legend:  IL  -  Points  in  ^  direction  JL  -  Points  in  rj  direction 


Rcc  -  Surface  mesh  Reynolds  number  Ad  -  Angular  spacing  (deg) 
c  -  Stretch  factor  at  surface 

Subscripts:  av  -  average  min  -  minimum 

max  -  maximum 

4.2  Results 

Computations  with  each  of  the  methods  are  compared  with  experimental  values  published  by 
Holden  et  al.  [46].  Results  from  the  MC  scheme  are  displayed  in  Fig.  3  where  the  unnormalized 
computed  surface  pressure  is  plotted  versus  9,  the  angle  measured  in  degrees  along  the  cylinder 
with  reference  to  the  stagnation  streamUne,  0  =  0.  Not  much  variation  with  mesh  resolution 
is  observed  in  the  quantitative  pressure  results.  The  experimental  measured  pressure  values 
display  significant  scatter  and  are  consistently  lower  than  the  computations.  A  similar  observation 
was  made  by  Prabhu  et  al..  The  computed  stagnation  point  pressure  compaies  very  well  with 
the  results  obtained  with  the  Rayleigh  supersonic  pitot  pressure  formula  [47],  4.20psia,  and  the 
computed  results  of  Prabhu  et  al.  [44].  Fig.  4  shows  the  unnormalized  computed  heat  transfer  (in 
Btujps)  versus  9.  Overall,  the  agreement  with  the  experimental  values  is  excellent  even  with 
the  coarsest  mesh  computed  which  has  a  very  high  surface  cell  Reynolds  number  (Table  4.1).  The 
stagnation  heat  transfer  agrees  very  well  with  the  Fay  and  Riddell  value  (as  quoted  by  Prabhu  et 
al.)  of  49.5  Btu/ ps.  It  should  be  noted  that  the  MC  method,  proposed  as  a  correction  to  the 
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Figure  3:  Surface  pressure  in  Mach  16  Bow  past  a  blunt  body  with  the  MC 


scheme 


original  Steger  Warming  (OSW)  method,  is  a  tremendous  improvement.  Results  from  the  OSW 
method  (not  shown  for  the  purposes  of  brevity)  indicate  an  overprediction  of  heat  transfer  of  500%, 
400%  and  250%,  respectively,  over  the  Fay  and  Riddell  value  on  Grids  1,  2  and  3  respectively. 

One  drawback  of  the  MC  algorithm  over  the  OSW  algorithm  is  the  tendency  to  develop  the 
anomalous  “carbuncle”  solutions  referred  to  earlier.  The  OSW  algorithm  does  not  display  any 
such  tendency  for  the  grids  computed.  Since  the  MC  algorithm  is  supposed  to  revert  to  the 
OSW  method  at  shocks,  this  is  more  appropriately  characterized  as  a  drawback  of  the  switching 
algorithm  utilized  in  the  present  research.  Indeed,  for  the  results  displayed  here,  the  carbuncle 
problem  was  eliminated  by  explicitly  applying  the  OSW  scheme  at  a  few  points  on  either  side 
of  the  shock,  i.e.,  overriding  the  pressure-based  switching  mechanism  described  earlier.  Fig.  5 
exhibits  the  Mach  number  variation  versus  distance  (measured  relative  to  the  cylinder  center) 
and  normalized  by  the  cylinder  radius  along  the  stagnation  streamline.  For  each  mesh,  the  shock 
is  captured  with  at  most  two  points  inside.  The  shock  standoff  distance,  defined  with  reference 
to  the  location  of  the  sonic  point,  is  uniform  with  mesh  resolution  (A/R  ~  0.4).  These  results 
are  in  close  agreement  with  the  results  of  Prabhu  et  al.  and  reflect  well  upon  the  shock  capture 
capability  of  the  OSW  method. 

For  the  sake  of  completeness,  we  show  the  computed  Mach,  pressure  and  temperature  contours 
in  Fig.  6.  With  the  strategy  of  explicit  application  of  the  OSW  scheme  in  the  vicinity  of  the  shock, 
no  evidence  of  the  carbuncle  is  visible. 

Fig.  7  shows  the  surface  pressure  prediction  with  the  van  Leer  method.  As  for  the  MC  results, 
the  pressure  prediction  does  not  change  much  with  grid  resolution  and,  in  fact  the  pressure  over 
the  entire  surface  of  the  cylinder  differs  only  by  an  insignificant  amount  from  that  predicted  by 
the  MC  method.  In  sharp  contrast,  grid  resolution  has  a  drastic  effect  upon  the  heat  transfer 
prediction  (Fig.  8).  Although  all  the  heat  transfer  profiles  are  smooth,  results  on  the  coarse  mesh 
overpredict  the  heat  transfer  by  a  large  amount  over  the  entire  surface  of  the  cylinder  (120%  over 
the  Fay  and  Riddell  value  at  the  stagnation  point).  The  results  improve  significantly  with  grid 
resolution  with  52%  overprediction  of  stagnation  heat  transfer  on  Grid  2  and  a  relatively  small  15% 
overprediction  on  Grid  3.  These  results  indicate  the  extreme  diffusiveness  of  the  van  Leer  scheme. 
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Figure  4:  Heat  transfer  in  Mach  16  flow  p<ist  a  blunt  body  with  the  MC  scheme 
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Figure  5;  Mach  number  along  stagnation  streamline  for  Mach  16  flow  past  a 
blunt  body  with  the  MC  scheme 
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a  result  previously  observed  by  other  researchers  on  other  geometries  and  flow  parameters  [22]. 
Despite  the  significant  inaccuracies  in  the  heat  transfer  prediction,  the  method  exhibits  excellent 
shock-capture  capability  (Fig.  9).  For  all  mesh  resolutions,  the  shock  is  not  only  captured  within 
two  zones  along  the  centerline  but  also,  the  shock  standoff  for  each  grid  is  within  one  grid  point. 
It  is  noted  that  none  of  these  calculations  with  van  Leer’s  scheme  tended  to  exhibit  any  tendency 
towards  the  development  of  the  anomalous  “carbuncle”  solution.  The  computed  Mach,  pressure 
and  temperature  contours  are  presented  in  Fig.  10. 

An  attempt  was  made  to  improve  the  predictions  of  heat  transfer  by  investigating  different 
choices  of  the  extrapolated  vector  W  (see  Section  2).  As  mentioned  earlier,  all  of  the  calculations 
described  in  this  work  utilize  W  =  {/?,  u,  u,p}  for  the  extrapolation  procedure  to  obtain  the 
L  and  R  states.  The  effect  of  utilizing  W  =  {/),pu,pi;,p}  and  the  set  of  conserved  variables 
^  was  investigated  for  the  van  Leer  scheme  on  Grid  2.  It  was  found  that  neither 

of  these  had  any  significant  influence  on  the  accuracy  of  the  final  result  (Fig.  11)  and,  in  fact,  the 
rate  of  convergence  was  found  to  be  adversely  affected. 

Turning  to  the  Roe  schema  (Fig.  12),  once  again  the  computation  of  the  surface  pressure 
displays  a  remarkable  insensitivity  to  grid  resolution.  This  is  perhaps  a  result  of  the  fact  that  the 
surface  pressure  is  a  mechanical  quantity,  dictated  mainly  by  the  shock  strength  and  shape  which 
in  turn  is  rather  accurately  predicted  by  all  the  inviscid  discretization  methods  at  each  level  of 
mesh  refinement.  In  this  regard,  the  fact  that  the  mesh  is  well  aligned  with  the  shock  wave  is 
an  important  factor.  For  Grids  2  and  3,  several  values  of  the  cutoff  parameter  6  were  utilized  for 
reasons  outlined  below.  The  effect  of  this  parameter  on  the  computed  surface  pressure  is  clearly 
negligible. 

Fig.  13  exhibits  the  surface  heat  transfer  comparison  with  experiment.  For  Grid  1,  the  value 
of  b  used  was  quite  small  (0.05).  No  tendency  toward  the  development  of  the  carbuncle  was 
observed  and  the  results  are  quite  accurate  in  similarity  with  the  predictions  with  the  MC  scheme. 
As  mentioned  previously,  Roe’s  scheme  exhibits  a  tendency  to  display  the  carbuncle  for  this 
configuration  as  observed  previously  by  Peery  et  al.  [31]  and  Liou  et  al.  [33].  In  the  present 
instance,  this  tendency  is  suppressed  by  increasing  the  value  of  the  constant  I  in  the  entropy 
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Figure  7:  Surface  pressure  in  Mach  16  flow  past  a  blunt  body  with  the  van  Leer 
scheme 
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Figure  9:  Mach  number  along  stagnation  streamline  for  Mach  16  flow  past  a 
blunt  body  with  the  van  Leer  scheme 
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(a)  Marh  contours  (b)  Pressure  contours  Temperature  contours 
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(b)  Surface  heat  transfer 

Figure  11:  Effect  of  extrapolation  of  different  variables  on  surface  quantities  with 
the  van  Leer  scheme 


29 


e 


Figure  12:  Surface  pressure  in  Mach  16  flow  past  a  blunt  body  with  the  Roe 
scheme 
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correction  formula  (Eqns.  2.26  and  2.27).  For  Grid  1,  we  believe  the  truncation  error  provided 
by  the  coarseness  of  the  mesh  is  possibly  sufficient.  For  Grid  2,  small  values  tended  to  give  the 
bulge  solution  clearly  visible  in  the  Mach  contours.  At  about  S  =  0.2,  the  visible  evidence  of  the 
bulge  in  the  Mach  contours  disappeared.  However,  a  small  dip  remained  in  the  surface  quantities 
of  interest  in  the  vicinity  of  the  stagnation  point.  Upon  increasing  S  to  0.3,  all  evidence  of  the 
carbuncle  disappeared.  As  shown  in  Fig.  13,  further  increase  to  0.4  did  not  have  any  significant 
effect  upon  the  flowfield.  Although  a  rigorous  examination  of  the  minimum  value  allowable  for 
an  anomaly  free  solution  was  not  made,  for  Grid  3,  the  value  of  S  was  increased  in  increments  of 
0.1  until  all  evidence  of  anomaly  in  contour  plots  of  the  Mach  number,  pressure  and  temperature 
disappeared.  A  minimum  value  of  0.5  was  required.  Even  then,  a  small  dip  persisted  in  the 
heat  transfer  computation  in  the  vicinity  of  the  stagnation  point  (Fig.  13).  Further  increase 
upto  S  =  0.8  did  not  remove  this  anomaly.  The  shock  standoff  distance  and  shock  capture 
capability  (Fig.  14)  are  similar  to  that  exhibited  by  the  van  Leer  scheme.  In  general,  the  shock 
is  captured  within  at  most  two  points  even  with  the  highest  value  of  S  utilized.  The  computed 
flowfield  exhibits  no  evidence  of  the  carbuncle  in  the  flow  contours  (Fig.  15).  Further,  the  standoff 
distance  (roughly  40%  of  the  radius)  is  quite  accurate  and  invariant  with  grid  resolution.  On  all 
schemes  considered,  this  reflects  very  well  upon  the  computed  density  ratio  across  the  bow  shock 
and  overall  low  conservation  errors. 

The  excellent  agreement  observed  for  heat  transfer  results  with  the  MC  and  Roe  scheme  may 
be  fortuitous.  Further  numerical  studies  are  necessary  at  different  Mach  and  Reynolds  number 
conditions  to  validate  the  conclusions  derived  in  the  previous  discussion. 
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Figure  13;  Heat  transfer  in  Mach  16  flow  past  a  blunt  body  with  the  Roe  scheme 
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Figure  14;  Mach  number  along  stagnation  streamline  for  Mach  16  flow  past  a 
blunt  body  with  the  Roe  scheme 
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5.  Corner  Flow 


5.1  Flow  Parameters  and  Xirid  Details 

The  parameters  for  this  flow,  chosen  to  simulate  the  experimental  evidence  of  Holden  et  al.  [34]), 
are: 

Moo  =  14.1 
Ramp  Angle  =  24° 

Too  =  130.8R 

Poo  =  0.21 

Reynolds  Number  =  72000per  foot 
^wall  = 

Two  types  of  grid  are  considered  for  this  study.  Grids  1,  2  and  3  are  are  simply  sheared  (Fig.  2(bl)) 
exponentially  stretched  on  either  side  of  the  interaction  and  also  in  the  surface  norma!  direc¬ 
tion.  Since  recent  research  efforts  have  highlighted  the  importance  of  resolving  the  leading  edge 
shock  [36],  Grid  4  utilizes  209  x  59  points  providing  better  resolution  in  the  vertical  direction 
at  the  leading  edge  by  suitable  redefinition  of  the  domain  at  the  leading  edge  (Fig.  2(b2)).  The 
horizontal  spacing  over  the  solid  surface  is  the  same  as  for  Grid  3.  The  characteristics  of  the  grids 
utilized  are  presented  in  Table  5.1. 

5.2  Results 

Fig.  16  (a),  (b)  and  (c)  display  the  comparison  of  pressure  coefficient  (defined  as  )  with 

the  MC,  van  Leer  and  Roe  schemes,  respectively.  The  abscissa  is  the  streamwise  distance  nor¬ 
malized  by  the  distance  from  the  leading  edge  to  the  corner  (L  =  1.44/f).  On  the  coarse  mesh, 
all  methods  predict  the  start  of  pressure  rise  at  an  Xf  L  value  of  1.0,  significantly  upstream  of  the 
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Table  5.1:  Corner  grid  details 


Grid 

IL  X  JL 

R^c\av 

Ax|a„ 

A 

49  X  14 

47.6 

0.02 

0.12 

0.04 

1.4 

B 

99  X  29 

21.7 

0.01 

0.06 

0.02 

1.2 

C 

199  X  59 

10.3 

0.005 

0.03 

0.01 

1.1 

D 

209  X  59 

5.9 

0.005 

0.03 

0.01 

1.1 

Legend:  IL  -  Points  in  ^  direction  JL  -  Points  in  t]  direction 

Rcc  -  Surface  mesh  Reynolds  number  Aa:  -  Streamwise  spacing  (ft) 
c  -  Stretch  factor  at  surface 

Subscripts:  av  -  average  min  -  minimum 

max  -  maximum 

experimental  value.  Further,  the  peak  pressure  location  as  well  as  magnitude  are  also  consider¬ 
ably  underpredicted.  For  Grid  2,  once  again,  no  significant  differences  between  the  methods  are 
observed  in  the  pressure  coefficent  prediction  though  some  improvement  is  observed  of  the  peak 
pressure  (underpredicted  roughly  22%)  as  well  as  the  initial  pressure  rise  over  that  obtained  on 
Grid  1.  In  contrast  to  Grids  1  and  2  where  not  much  difference  is  observed  between  the  schemes, 
the  pressure  coefficient  on  Grid  3  displays  significant  differences  between  the  methods.  The  initial 
pressure  rise  occuring  at  about  xj L  ~  0.6  is  best  reproduced  by  the  MC  method  while  for  the 
Roe  and  van  Leer  schemes,  this  rise  occurs  slightly  downstream  {xj L  ~  0.75).  The  location  of 
the  peak  pressure  is  also  best  captured  with  the  MC  method  although  the  magnitude  is  overpre¬ 
dicted  about  12%.  On  the  other  hand,  the  location  of  the  computed  peak  with  the  Roe  and  van 
Leer  schemes  is  slightly  upstream  of  the  experimental  value  though  the  magnitude  is  correctly 
predicted.  There  is  not  much  difference  in  the  computed  pressure  coefficients  with  Grid  4  over 
Grid  3  for  both  the  MC  and  Roe  schemes.  However,  the  peak  Cp  with  van  Leer’s  scheme  drops 
somewhat  (3%).  All  schemes  predict  the  theoretical  inviscid  pressure  rise  beyond  xj L  ~  2  on 


36 


every  grid. 

The  heat  transfer  results  are  displayed  in  Fig.  17  (a),  (b)  and  (c)  for  the  MC,  van  Leer  and 
Roe  schemes,  respectively.  The  ordinate  in  this  figure  is  the  heat  transfer  coefficient  defined  as 
^  (in)  IPoo'Uqo  {Ho  —  H^)  where  H  is  the  total  enthalpy  and  the  subscript  w  denotes  evaluation 
at  the  wall.  The  inaccuracies  on  the  coarsest  mesh  may  be  predominantly  attributed  to  inadequate 
grid  resolution  which  leads  to  significant  truncation  error.  For  Grid  2,  both  the  MC  method 
and  Roe’s  method  display  higher  peak  heat  transfer  rates  while  van  Leer’s  method  exhibits  a 
lower  peak,  each  method  approaching  closer  to  experimental  values  relative  to  Grid  1.  On  the 
finest  mesh,  once  again,  the  MC  method  predicts  the  initial  drop  (corresponding  actually  to  the 
point  of  separation  as  evident  feter)  at  the  correct  location  while  the  Roe  and  van  Leer  schemes 
underpredict  upstream  influence.  On  the  other  hand,  the  peak  value  is  best  predicted  with  Roe’s 
scheme  (15%  overprediction)  with  the  MC  method  overpredicting  the  heat  transfer  coefficient 
the  most  (30%  overprediction).  In  a  manner  similar  to  that  observed  previously  for  the  pressure 
coefficient,  there  is  no  improvement  with  Grid  4  over  Grid  3.  Note  that  with  Roe’s  scheme,  a 
small  dip  is  observed  at  Xf  L  ~  1.25  on  Grids  1  and  2,  the  cause  of  which  is  presently  unknown. 
However,  with  grid  refinement,  this  anomaly  is  eliminated. 

Fig.  18  compares  the  prediction  of  the  skin  friction  coefficient,  defined  as  2/iu,  fir) 

The  size  of  the  separated  reverse  flow  region  is  indicated  by  the  region  of  negative  C/.  This 
increases  with  grid  refinement  for  all  methods.  On  Grids  1  and  2,  once  again  the  MC  and  Roe 
schemes  yield  similar  results  and  van  Leer’s  method  overpredicting  significantly  the  peak  skin 
friction.  All  methods  predict  almost  negligible  separation  on  Grid  1  indicating  a  lack  of  adequate 
resolution  of  the  viscous  terms.  For  the  denser  Grids  (3  and  4),  the  predictions  with  van  Leer 
algorithm  rapidly  approach  those  with  the  other  two  schemes,  a  trend  also  observed  earlier  with 
the  heat  transfer.  The  MC  method  predicts  the  most  accurate  start  of  the  separation  region  in 
comparison  with  experimental  values.  All  three  methods,  however,  fail  to  correctly  locate  the 
trough  in  skin  friction  coefficient,  most  likely  due  to  reattachment,  at  X/L  ~  1.25. 

Overall,  the  effect  of  grid  resolution  on  all  methods  is  to  provide  better  comparison  between 
theory  and  experiment  as  anticipated.  For  the  sake  of  completeness,  the  computed  Mach  contours 
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on  Grid  4  with  each  method  are  presented  in  Fig.  19.  The  overall  pattern  is  the  same  for  all 
methods  and  similar  to  that  observed  on  Grid  3  (not  shown).  The  leading  edge  shock  is  captured 
accurately  but  tends  to  get  smeared  near  the  point  where  it  intersects  the  corner  shock.  The 
resolution  of  the  flow  at  this  point  of  intersection  may  be  critical  for  the  purposes  of  accuracy. 
It  is  interesting  to  note  that  although  van  Leer’s  algorithm  performs  poorly  on  the  coarse  mesh, 
significant  improvement  is  obtained  with  grid  refinement  sepecially  in  the  prediction  of  surface 
quantities.  The  salient  difference  between  Meshes  3  and  4  is  in  the  vertical  spacing.  However, 
the  horizontal  spacing  is  identical  except  upstream  of  the  leading  edge  of  the  plate  where  Grid  4 
has  10  uniformly  spaced  points.  Clearly  the  issue  of  horizontal  spacing  has  not  been  resolved  and 
possibly  even  better  comparison  may  be  obtained  with  each  method  if  the  regions  around  the 
points  of  separation  and  reattachment  are  further  resolved. 

It  may  be  mentioned  that  the  present  results  on  Grids  3  and  4  agree  in  an  overall  sense  with  the 
computations  of  other  researchers.  Although  precise  comparisons  are  not  possible  due  to  different 
implementations  and  mesh  details,  some  differences  in  trends  may  be  outlined.  The  CFL3D  code 
(utilizing  Roe’s  scheme  with  third  order  approximation  and  the  thin  layer  equations  for  the  viscous 
terms)  has  been  utilized  for  the  same  configuration  by  Rudy  et  al.  [36]  as  well  as  by  Rizzetta 
et  al.  [35].  The  results  of  the  former  indicate  surprisingly  good  agreement  with  experimental 
data  on  a  relatively  coarse  mesh  (51  X  51)  though  the  predictions  deteriorate  somewhat  with 
higher  grid  re.solution.  Their  grid  independent  results  (obtained  on  a  101  X  101  mesh)  exhibits 
a  larger  separation  region  than  that  observed  experimentally.  However,  they  note  that  their 
time-accurate  calculations  require  about  12ms  for  steady  state  to  be  achieved.  Sms  more  than 
observed  experimentally  and  2ms  more  than  the  total  experimental  time.  Indeed,  their  best 
comparison  with  experimental  results  is  obtained  during  the  transients  (between  2  and  3ms)  in 
time-accurate  computations.  In  comparison,  with  the  convergence  criteria  outlined  earlier,  the 
results  displayed  on  the  Mesh  4  represent  a  minimum  flow  development  time  of  9.5ms  for  van 
Leer's  scheme  the  relevant  numbers  for  the  MC  scheme  and  Roe’s  scheme  are  19ms  and  17ms, 
respectively.  One  difference  in  trends  is  that  the  present  methods  all  display  pressure  and  heat 
transfer  peaks  upstream  of  the  experimentally  observed  location  in  contrast  to  those  obtained 
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Corner  flow  (Grid  4)  -  MC 


Corner  flow  (Grid  4)  -  vL 


Corner  flow  (Grid  4)  -  Roe 


Figure  19;  Mach  contours  for  compression  corner  flow  on  Grid  4 


by  Rudy  et  al.  for  whom  the  peaks  lie  downstream  of  the  experimental  locations.  The  results 
of  Rizzetta  et  al.  with  the  CFL3D  code  indicate  a  higher  sensitivity  to  vertical  grid  resolution 
than  observed  in  the  present  research.  Results  with  the  Beam- Warming  algorithm,  the  VAPOR 
code  (based  upon  the  OSW  algorithm)  and  MacCormack  explicit  predictor  corrector  algorithm 
were  also  presented  by  Rizzetta  et  al..  On  a  grid  identical  to  that  denoted  Grid  3  in  the  present 
research,  the  Beam- Warming  and  OSW  algorithms  exhibit  significant  underprediction  in  peak 
heat  transfer  while  the  MacCormack  predictor  corrector  algorithm  performs  reasonably  well.  All 
the  computations  mentioned  above,  as  also  those  reported  by  Thareja  et  al.  [48]  with  the  adaptive 
LAURA  algorithm  [48]  fail  to  predict  the  dip  in  skin  friction  coefficient  at  about  Xj L  ~  1.25. 
The  location  of  separation  and  reattachment  points  and  extent  of  separation  are  compared  in 
Table  5.2  for  several  research  efforts. 
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Table  5.2:  Separation  region  details  for  corner  flow 


Solution 

Separation 

Reattachment 

Total 

Algorithm 

Point 

Point 

Extent 

Experiment 

0.5475-0.7025 

1.1278-1.3748 

0.5703-0.8273 

CFL3D  [48] 

0.5318 

1.3431 

0.8113 

Thareja  et  al.  [48] 

0.5588 

1.3493 

0.8113 

A 

0.9813 

1.0242 

0.0429 

MC 

B 

0.8929 

1.1100 

0.2171 

C 

0.6516 

1.2748 

0.6232 

D 

0.6596 

1.2688 

0.6092 

■ 

A 

— 

— 

No  sep.  detected 

D 

B 

0.9274 

1.1244 

0.1970 

■ 

C 

0.7703 

1.2150 

0.4446 

■ 

D 

0.7733 

1.2166 

0.4432 

A 

0.9871 

1.0173 

0.0302 

Roe 

B 

0.9033 

1.1226 

0.2194 

C 

0.7449 

1.2367 

0.4917 

D 

0.7330 

1.2381 

0.5051 
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6.  Conclusions 


A  detailed  examination  of  the  performance  of  several  flux-split  algorithms  is  attempted  on  two 
problems  representative  of  complex  flowfields.  Results  are  compared  with  experimentally  observed 
data.  Based  upon  available  results,  the  following  conclusions  are  made: 

•  The  method  of  MacCormack  and  Candler  performs  accurately  for  both  the  blunt  body  flow 
as  well  as  the  ramp  flow  insofar  as  both  surface  pressure  and  heat  transfer  are  concerned. 
The  overall  accuracy  is  comparable  to  Roe’s  flux-difference  split  scheme  and,  in  the  present 
numerical  framework,  is  more  robust  -  permits  code  operation  at  higher  CFL  numbers 
leading  to  more  rapid  convergence. 

•  van  Leer’s  scheme:  This  scheme  displays  a  tendency  to  consistently  overpredict  heat  transfer 
indicating  perhaps  excess  diffusion  in  the  boundary  layer.  However,  with  grid  refinement, 
the  accuracy  of  this  method  rapidly  approaches  those  of  the  other  two. 

•  Roe’s  flux  difference  scheme:  This  scheme  exhibits  a  tendency  toward  the  development  of 
anomalous  “carbuncle”  solutions  for  blunt  body  flows.  This  tendency  may  be  suppressed 
by  appropriate  increase  in  the  eigenvalue  cutoff  required  to  enforce  the  entropy  condition. 
For  the  corner  flow,  the  performance  of  this  scheme  is  comparable  to  the  MC  scheme. 

Although  the  MC  and  Roe  schemes  exhibit  excellent  comparison  with  experimental  data,  this 
result  is  representative  only  at  the  flow  parameters  considered  and  may  not  carry  over  to  other 
conditions.  Further  computations  at  higher  Mach  and  Reynolds  numbers  are  necessary.  For 
corner  flows,  the  present  results  indicate  the  need  for  further  computations  with  higher  streamwise 
resolution. 
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A 

B 

Btu 

c 

C 

C, 

CFL 

e 

ei 

f 

F, F,G,G 

G. R. 
IL,JL 
J 

L 

M 

MC 

n 

Hit 

osw 

p 

psi 

Pr 

PR 


Nomenclature 

flux  Jacobian  of  F 

flux  Jacobian  of  G 

British  Thermal  Units 

stretch  factor;  local  speed  of  sound 

component  of  Q 

specific  heat  at  constant  volume 
Courant- Friedrich- Levy  number 
total  energy 
internal  energy 
feet 

flux  vectors 

global  residual 

points  in  ^  and  rj  direction 

Jacobian,  inverse  ceU  volume 

distance  from  leading  edge  to  corner 

Mach  number 

MacCormack  and  Candler  scheme 
normal  direction 

number  of  iterations  to  double  CFL 
Original  Steger  Warming  scheme 
pressure 

pounds  per  square  inch 
Prandtl  number 
pressure  function 
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Q 

matrix  diagonalizing  A 

R 

gas  constant;  Rankine 

Re 

Reynolds  number 

RMSP 

root  mean  square  surface  pressure 

RMSQ 

root  mean  square  surface  heat  transfer 

S 

Jacobian  of  V  with  respect  to  U 

t 

time 

T 

temperature 

u 

Cartesian  velocity  in  x  direction 

U,U 

solution  vector 

V 

Cartesian  velocity  in  y  direction 

V 

Cartesian  velocity  vector 

vL 

van  Leer  scheme 

V 

vector  of  primitive  variables 

W 

vector  of  extrapolated  variables 

Cartesian  coordinates 

6 

cutoff  value;  any  change 

S 

cutoff  parameter 

A 

shock  standoff  distance 

V 

gradient  operator 

7 

ratio  of  specific  heats 

A 

eigenvalue 

d 

partial  derivative  operator 
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p 

9 

Xi 


density 

angle 

transformed  coordinates 
pressure  switch 


Subscripts 

j  +  h 


V 

w 

oo 


interface  between  j  and  j  +  1 

viscous 

wall 

freestream  conditions 


Superscripts 


L 

R 


+  » ~ 


state  at  left  of  interface 
state  at  right  of  interface 
positive  and  negative  components 
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